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A new code to study structures in collisionally active, perturbed 
debris discs. Application to binaries. 
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ABSTRACT 

Context. Debris discs are traditionally studied using two distinct types of numerical models: statistical particle-in-a-box codes to 
study their collisional and size distribution evolution, and dynamical N-body models to study their spatial structure. The absence 
of collisions from N-body codes is in particular a major shortcoming, as collisional processes are expected to significantly alter the 
results obtained from pure N-body runs. 

Aims. We present a new numerical model, to study the spatial structure of perturbed debris discs at dynamical and collisional steady- 
state. We focus on the competing eff'ects between gravitational perturbations by a massive body (planet or star), collisional production 
of small grains, and radiation pressure placing these grains in possibly dynamically unstable regions. 

Methods. We consider a disc of parent bodies at dynamical steady-state, from which small radiation-pressure-aff'ected grains are 
released in a series of runs, each corresponding to a diff'erent orbital position of the perturber, where particles are assigned a collisional 
destruction probability. These collisional runs produce successive position maps that are then recombined, following a complex 
procedure, to produce surface density profiles for each orbital position of the perturbing body. 

Results. We apply our code to the case of a circumprimary disc in a binary. We find pronounced structures inside and outside the 
dynamical stability regions. For low Cb, the disc's structure is time varying, with spiral arms in the dynamically "forbidden" region 
precessing with the companion star. For high ^5, the disc is strongly asymmetric but time invariant, with a pronounced density drop 
in the binary's periastron direction. 
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1. Introduction 

Debris discs have been detected around a significant fraction (be- 
tween ~ 5 and ~ 3 0%) of main sequence stars of all spectral 
types (e.g. lSu etall [20Q6: Trilling etal. 2008). They consist of 
collisionally interacting solid bodies spanning a wide size range, 
from -1-100 kilometre- sized parent bodies, probably leftovers 
from the planet formation process, down to micron-sized debris. 
As they make up most of the geometrical cross section, only 
the smallest, < 1 cm, particles of this collisional cascade are de- 

. tectable, in t hermal emission or scattered light (see review by 

■ IWvattll2008h . 

] For most discs that have been resolved, pronounced 
structures have been observed, such as t wo-sided asymme- 
tries, spirals, warps, clumps or rings (e.g [Kalas et al.l 120051: 

; iGolimovski et al.ll2006l:rSchneider et S]l2009l) . Such features in- 
dicate complex dynamical environments, and numerous studies 
have been aimed at explaining their origin (e.g. Krivov 2010). 
Most of these studies are based on A/^-body numerical codes that 
follow the dynamical evolution of the system, exploring, de- 
pending on each system's specificities, diflTerent scenarios: per- 
turbing influence of (often undetected) planets or stellar compan- 
ions, violent transient events, or interaction with gas remnants. 
In such a collisionless A/^-body approach, the disc is sampled by 
a population of Nnum test particles whose trajectories can be pre- 
cisely tracked. The forces that control their evolution being: the 
central star's gravity, the gravitational pull of planets and other 
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stars, radiation pressure, Poynting-Robertson (PR) drag, stellar 
wind drag. The precision of the N-ho&y method allows one to 
study fine efifects such as orbital resonances or transient spa- 
tial inhomogeneities. However, this comes at a price, which is 
that collisions are neglected. This leads to several shortcomings, 
which can be more or less worrying depending on the issues that 
are being addressed. These shortcomings can be the following: 

- l)No size distribution evolution. Since collisions control the 
particle size distribution and its evolution (by merger, ero- 
sion or shattering), collisionless N-body codes cannot handle 
this important issue. In particular, they cannot assess to what 
extent the specific dynamics of a given system can aflTect the 
physical evolution of its population, since the number, size 
distribution and velocity spread of collisional fragments di- 
rectly depend on impact velocities, and thus how size distri- 
butions can vary as a function of spatial location. 

-2) No feedback of collisions on the dynamics. Impacts, be 
it accreting, bouncing or fragmenting ones, necessarily alter 
the orbits of colliding objects, because angular momentum is 
redistributed and kinetic energy is lost to heat. These efiTects 
can significantly change the results of pure N-body simula- 
tions by damping high eccentricities induced by a planetary 
or stellar perturbations, or by ejecting collisional fragments 
far from their progenitors. 

-3) Particle lifetimes / Spatial structures. Dynamical pro- 
cesses have characteristic timescales, which can be rela- 
tively long, as for instance for mean motion resonances 
or secular eflTects. Depending on the disc's density, colli- 
sion timescales could be shorter than the dynamical ones, 
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thus interfering with, or even preventing the development 
of the corresponding spatial structures. This problem gets 
even more complex because collisional timescales are in 
most cases size-dependent. For the smallest grains, of spe- 
cific interest because they usually dominate the disc luminos- 
ity (iThebault & Augereaull2QQ7h , the dynamical timescales 
also become size-dependent because of the eff'ect of radia- 
tion pressure. These competing timescale eff'ects, in addition 
to aff'ecting the development of spatial structures, can have 
direct consequences on the size distribution, which can, at 
some locations in the disc, strongly depart from any classi- 
cal equilibrium power-law. 

Incorporating collisions into an N-body scheme is extremely 
difficult in the context of debris discs, where typical impact ve- 
locities are high and lead to destructive collisions producing a 
vast amount of small fragments. Following all these fragments 
would lead to an exponential increase in the number of particles 
that would soon become unmanageable. 

The collisional evolution of debris discs is usually investi- 
gated separately with a different type of model, based on sta- 
tistical particle-in-a-box codes, where the solid body popula- 
tion is distributed into logarithmically- spaced size bins, whose 
number density is followed using detailed prescriptions for col- 
lision outcomes (fragmentation, cratering, accretion). The price 
to pay here is a poor spatial resolution and a very limited mod- 
elling of the dynamics: collision rates are estimated using spa- 
tially averaged estimates of orbital paramet ers that are fixed or 
follow very simplified evolution laws (e.g. iKrivov et"al 
iThebault & Augereau ' 2007). A pure statistical approach thus 
cannot give any precise information on the creation and evolu- 
tion of spatial structures. For this, an N-body scheme is unavoid- 
able. 

1.1. N-body and collisions: Previous studies 

Given the sheer difficulty of the task, there have only been a 
few attempts at incorporating fragmenting collisions into N- 
body codes. The conceptually most natural, but in practice most 
difficult way to do this is by a "brute force" method that fol- 
lows the fragments produced after each impact. The pio neering 
study of that method is that of iBeauge & AarsethI (Il990h . whose 
code followed 4 fragments after each shattering collision, but 
had to be restricted to only 200 initial bodies, and was limited 
to relatively short timescal es before reaching a crit i cal nu mber 
of particles. More recently. ILeinhardt & RichardsonI (l2005l) con- 
sidered the breakup of rubble piles of gravitationally bound hard 
spheres, but did not get below the size of the initially defined 
"hard" sphere units. 

Another solution is to use a mix of the N-body and the statis- 
tical approaches. To our knowledg e, the on ly published attempt 
in this direction is that of Grigorie va et al.l |2007), using a pop- 
ulation of "super particles" (SP), whose orbits are determinis- 
tically followed, each representing a cloud of monosize grains. 
When SPs' paths cross, collisions are considered between their 
respective grain populations in a statistical way, and new SPs 
are created to account for the fragments. However, this model 
was restricted to very short timescales of a few orbital periods. 
This combined N-body/statistical approach is thus far from be- 
ing able to address long term phenomenon, but it is probably the 
most promising one in terms of its potentialities in solving all 3 
main problems listed in the previous section. 

As of today, t he be st available model is probably that of 
IStark & Kuch'i^ (l2009h . whose main ambition is to address 



shortcoming number 3 Q regarding the competing effects of dy- 
namical and collisional lifetimes on the development of spatial 
structures. Its principle is to release a population of collision- 
less test particles and record, at regularly spaced time intervals, 
their position to construct a stream of positions and velocities 
for each particle. All streams are then recombined to create a 
surface density map on which the same particles are released 
once again and removed following collision probabilities derived 
from the density map. A new map of streams is then created and 
the process is iterated until it converges. This pioneering model 
has gi ven impressive results for the specific case of the Kuiper 
Belt ( Kuchner & Stark"2010). In its current version, this code is 
designed for the specific case of one perturber on a circular orbit. 
Although the case of multi-perturbers and of an eccentric per- 
turber could in principle be implemented, these issues have not 
been addressed yet and the code might not be able to handle fast 
evolving spatial structures (Stark, private communication). In the 
case of an eccentric perturber, this code requires that the time be- 
tween two consecutive records must equal a multiple of the or- 
bital timescale, which must be less than the collision timescale. 
As a result, this code may be poor at modelling highly collisional 
discs with distant perturbers on eccentric orbits, such as the case 
of a circumprimary debris disc in a binary. 

We present here a new code, also aimed at studying how 
collisional lifetimes affect the development of perturbed spatial 
structures (issue number 3), which can be used for the generic 
case of a collisional disc perturbed by one massive body, planet 
or stellar companion, having any possible orbit, circular or ec- 
centric, internal to, embedded in, or external to the disc. This 
code can handle size-dependent effects, in particular for small 
grains placed on high-eccentricity orbits by radiation pressure. 
The only required assumption is that a dynamical and collisional 
steady state has been reached in the system. We illustrate this 
new model with the specific case of a circumprimary debris disc 
perturbed by an external stellar companion. 

2. Model 

We assume that the studied system has reached a steady state, 
both dynamically and collisionally. We consider the following 
forces: gravity and radiation pressure from the central star, grav- 
itational pull from a perturbing planet and/or stellar companion. 
Collisions are assumed to produce bodies following a size dis- 
tribution in dN oc s^ds. We divide particles into two categories: 
parent bodies (PB), i.e., bodies that are large enough not to be 
affected by radiation pressure, and small fragments steadily pro- 
duced by collisions from these parent bodies and placed on ec- 
centric or unbound orbits by radiation pressure. The magnitude 
of the radiation pressure force is quantified, as usual, by the pa- 
rameter 13, corresponding to the size-dependent ratio of this force 
to that of star's gravity. The dynamical evolution of the parent 
bodies and the smaller grains is computed using an adaptive- step 
4^^ order Runge-Kutta integrator. 

For the sake of simplicity, we consider a reference case of 
particles orbiting a central star and one external perturber on an 
eccentric orbit. Our numerical procedure is divided into three 
steps that can be schematically described as follows. 

^ Although not implemented yet, some degree of fragmentati on is 
in principle possible to incorporate in this algoritm ( Stark & Kuc hned 
2009), but probably not enough to fully address points number 1& 2, 
especially the way impact velocities locally affect size distributions and 
velocity spreads by controlling the production of clouds of small frag- 
ments 
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Fig. 1. Schematic presentation of the numerical model. Usav = 10 "collisional" runs are performed, each taking as a starting point 
one of the Usav parent body (PB) disc configurations, corresponding to Usav diff'erent positions of the companion star on its orbit 
separated by a constant time interval A^^^y = torbl^sav^ obtained at the end (i.e., steady state) of the parent body run. Nnum = 2x10^ 
particles are released from the PB's positions, following a size distribution dN oc s^ds down to the radiation pressure blow out 
size SRp. Particles positions are recorded at time intervals separated by /^tsav Each simulation is run until all particles have been 
removed, either by dynamical encounter with the companion star or by collisions. The collected data is then used to reconstruct the 
dust distribution at dynamical and collisional steady state by means of the following procedure (indicated by the red arrows on the 
graph): the dust distribution at a given time to, corresponding to a given position of the companion star, is the combination of the 
dust released at ^ for the run, plus the dust particles released at to - Atsav for the - 1 run that have not been removed at to, plus 
the dust released at to - 2 Atsav for the - 2 run that has not been removed at to, etc. The procedure is iterated until we reach the 
j final record (particles released at - jfinai^tsav)^ for which no particle has survived until ^o- Given that the binary orbit is divided 
into Hsav = 10 positions, all - 10 x j runs correspond to the one. (The diff'erent disc profiles displayed in the figure are not 
simulation results but only illustrative sketches). 
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2.1. Parent Body Runs 

A first run is performed with only the parent bodies and thus 
only the central and perturbing gravitational forces. This run is 
stopped once a dynamical steady state has been reached, i.e., 
when no changes are observed between two consecutive pas- 
sages of the perturber at the same phase. At this steady state, the 
position of the parent body disc is then recorded for a series of 
Hsav diff'erent positions of the perturber on its orbit, separated by 
a fixed time interval /^tsav = tporbl^sav (where tporb is the orbital 
period of the perturber). 

2.2. CoHisional runs 

A series of "collisional" runs are then performed, each taking 
as an initial condition the positions of the parent bodies for one 
of the Hsav phases of the perturber, as indexed by the number 
i(p, with I < i(p < Hsav, from the sample of recorded PB con- 
figurations. For each collisional run, at ^ = 0, Nnum particles 
are releas ed following a s ize distribution in dN oc s^ds, with 
q = -3.5 (iDohnanyill 19691) . The positions of these particles are 
then recorded, at regularly space time intervals Atsav correspond- 
ing to the same sample of orbital phases of the perturber as in the 
PB run. At each integration timestep, each of the released par- 
ticles is assigned a collisional destruction probability focoii de- 
pending on its size s, velocity and the local geometrical optical 
depth Tr,e of the system: 

focoll = — TT7 — (1) 

where dt is the time increment of the code, torb the orbital period 
of a /5 = body at this location of the disc, is a reference 
size, dV and ^Vq are the departures from the local Keplerian 
velocity Vkcp of the conside red particle and for a parent body 
with 13 = respectively (see iThebault et al.ll201Ql for more de- 
tails on the procedure). The optical depth values T(r,e) are de- 
rived from the risav parent body runs, under the assumption that 
collisions occur mainly in the parent body region. This assump- 
tion is in agreement with previous debris disc niodell ing results 
(e.g. iKrivov et al.1 12Q06I: iThebault & Augereau I [200 7) and has 
been adopted i n all studies investigating the outer regions of de- 
bris discs (e.g lSmibbe_&Chiang 2006; Thebaul t & Wu 2008). 
In practice, we record the (steady-state) final positions of all par- 
ent body runs and transform them into (r, 6) maps of the vertical 
optical depth. These maps only give relative values and have to 
be normalized by (r), the average normal optical depth assumed 
for the PB disc. The collisional activity can be thus tuned in by 
setting the initial value for (r). Note that this procedure requires 
that dt is smaller than the collision time for the grain to properly 
resolve the collisional destruction probability, a criteria that is 
always met for the study of perturbed debris discs, where dt is a 
fraction of the orbital period and tcoii is always > tgrb- 

The collisional runs are stopped when all particles have been 
eliminated, either by collisions or by dynamical ejection after 
an encounter with the perturber. At the end of this step, to each 
initial phase of the perturber is assigned a series of surface 
density maps cr\i(p, 0), cr'(/0, 1), .., cr\i(p, j), recorded at the same 
regularly spaced time intervals Atsav, following the fate of all 
particles released when the perturber was at the position. Note 

^ In the case, not studied here, of an interior perturber, this method 
might not sample the whole orbit of a parent body. But this is not a 
problem as long as the number of bodies is large enough to populate all 
of given body's orbit longitudes. 



that all snapshots o-\i(p, j), cr^i^p, j -h risav), cr'i^, j -h 2nsav), etc., 
correspond to successive passages of the perturber at the same 
orbital phase (each perturber orbit being divided into risav posi- 
tions). 

2.3. Recombining 

In the final stage, we use all data collected in step 2 to recon- 
struct the dust distribution, at steady state, for each possible 
orbital phase of the perturber. The principle is that, at a given 
time ti^ corresponding to the perturber' s position Z^, the total dust 
population regroups grains that have just been produced (for the 
present position of the perturber), as well as survivor grains 
that have been produced earlier (when the perturber was at a dif- 
ferent orbital position), and that have not been collisionally de- 
stroyed or dynamically ejected yet. The procedure is the follow- 
ing: we start with the most recent dust particles, produced at ti^, 
whose spatial distribution is given by the saved record cr^i^p, 0). 
We then add the previous generation, produced at ti^- Atsav when 
the perturber was at angular position index - 1 , whose spatial 
distribution at time ti^ (i.e., at time Atsav after their release) is 
given by the saved record cr^i^p - 1,1). This procedure is then 
iterated, working our way back in time and piling up all the sur- 
viving grains from the successive records cr\i(p-j, j), to produce 
the total geometrical optical depth at time ti^ 

j^jmax 

cr(i<p) = ^ cr\i^- jj) (2) 

The index j^ax corresponds to the most ancient record, for 
which all initially released particles have been eliminated after 
the time jmax x ^tsav separating their release from the present 
time. Note that all the snapshots cr\i(p - 7, j), o-\i(p - j- risav, j + 
^sav), o-\i(f) - j - '^^sav, j + '^^sav), ^tc, are taken from the same 
source collisional run, with the perturber position at release be- 
ing (i(p-j), but separated by an integer number of orbital periods. 
In other words, the cr\i(p-j-k. risav, j^k.risav) record corresponds 
in practice to the o-\i(p - j, j -h k.risav) one. 

An illustrated summary of the main steps of our numerical 
procedure can be found in Fig[T] 

2.4. Tests 

As there is no analytical solution to the full problem of a col- 
lisionally active disc perturb ed by a massive body, w e follow 
an approach similar to that of IStark & Kuchnej (l2009l) and first 
check the correctness of our algoritm's solution for a simplified 
cas e with no companion. F or this case, Strubbe & Chiang (2006,) 
and lThebault & Wul (l2008h have shown that, for a collisional de- 
bris ring whose parent bodies distribution has a sharp outer edge 
rout, the radial profile of the surface density beyond rout should 
tend towards an aymptotic slope in r"^-^ (because of the small 
grains produced within the ring and pushed there on eccentric 
orbits). No te that we do not try to re pro duce the numerical sim - 
ulations of IStrubbe & C hiang (' 2006h or IThebault &Wul (l2008h . 
which did address slighly diff'erent issues, but the robust conclu- 
sions from their analytical derivations showing that a collisional 
ring always tends towards a profile in -1.5 outside the collisional 
active region. We thus run a test simulation with no companion 
and a parent body ring having a sharp outer edge at rout = 1 (ar- 
bitrary unit). As is clearly seen in Fig 12^, we do obtain a profile 
that tends towards the standard slope in -1.5 beyond r ~ 1.5. 
Note that this test is not a validation of the whole collisional 
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Fig. 2. Test runs: Left Panel. Azimutally averaged radial profile of the normalized surface density for a fiducial case with no 
companion. The expected theoretical asymptotic behaviour in r"^-^ is plotted as a reference. The initial extension of the parent body 
ring is set to r = 1 Right Panel: Azimutally averaged radial profile of the normalized surface density for a companion of eccentricity 
ei) = 0.5 as compared to the results of The bault et al.l (12010) with a diff'erent code. For this with-companion run, disntances have 
been normalized to the theoretical limit for orbital stability rent (see Sec B.lb . 



procedure, as the -1.5 slope is a result that does only weakly 
depend on the collisional rate within the PB ring. However, it 
validates the recombination procedure, which is rather complex, 
piling up maps after maps coming from diff'erent runs and at 
diff'erent times, none of these individual maps having a radial 
profile tending towards -1.5. 

For the case with a massive exterior companion, we perform 
an additional test by com paring our results to the ones obtained 
bv lThebaulte"taD(l2010h who, using a completely diff'erent algo- 
rtim with only 1-D resolution, estimated the azimutally averaged 
radial profile of the surface density. We take as a reference the 
case with a M2 = 0.25 Mi companion on a = 0.5 eccen- 
tric orbit and a parent body ring of optical depth r = 2 x 10"^ 
extending out to the orbital stability limit rcru (see Sec B.lb . We 
then azimutically average the spatial distribution and display it in 
Fig|2l As can be seen in Fig|2t) the obtained averaged radial pro- 
file is a very good match to the one obtained by ^Thebault et aP 
(I2OIOI) with their diff'erent method (see Fig. 2 of that paper), i.e. 
a steeper profile than in the previous no-perturber case because 
of the steady removal of high-jS particles by the companion. This 
second test is a robust validation of the collisional and dynamical 
prescriptions, as the shape of the radial profile directly depends 
on the balance between the collisional activity in the PB disc and 
the dynamical ejection of particles perturbed by the companion. 

3. A case study: debris disc in a binary 

To illustrate the potential of our code, we consider a case study 
of a circumprimary debris disc in a binary system, of particular 
interest in light of the recent surge of research regarding planet 
formation in binaries (for a recent review, see .Thebault 11201 lb . 
This case ha s been investigated in several recent works, both ob- 
serv ational (Trilling et al. 2007; Plavchan et al.""2009|; iDuchend 
I2OIO) and theoretical (Thebault et al. 2010). The main issue 
investigated in these studies was the extent of circumprimary 
discs, in particular if the companion star can induce a trunca- 



tion that can be detectable when looking at infra-red excesses 
or at the radial profile of the resolved disc. The numerical study 
of Thebault et al. ( 2010) showed that the coupling of collisional 
activity and radiation pressure plays a crucial role, steadily plac- 
ing small dust grains in regions that are in principle dynamically 
unstable. Since these grains dominate the flux at all wavelengths 
up to mid-IR, debris discs can thus appear to extend far beyond 
the theoretical radial distance rent for orbit al stability around 
the primary. However, lThebaultetal.1 (l2010l) used a collisional 
code with only 1-D spatial resolution (all azimutal information 
being lost in phase averaging), and were thus unable to study 
how binarity afl'ects the shape of circumstellar discs. This is- 
sue is crucial, as the presence of a companion star has been 
invoked as a po ssible explanation for several systems' aspe ct, 
such as HR4796 ("Ausereau et al." 1999; "Schneider et al.'"2009') or 
HD141569 feugereau & Papaloizou 2004: Ouillen et al. 2005.: 
lArdila et al.ll2Q05D . 

3.1. Setup 

We consider a binary of separation ag, eccentricity eB and mass 
ratio yu. We make the not- so-unreasonable assumption that the 
disc of large parent bodies has been shaped and truncated by the 
companion star. To ensure this, we allow the initial disc to extend 
slightly beyond the empirical ( 1-D) outer lim it rent for orbital 
stability derived by Holman & WiegertI (Il999h . and let the code 
naturally remove all unstable particles. 

To allow easy comparison between difl'erent cases, all dis- 
tances are normalized so that rent = 1 . In order to avoid the huge 
computational cost of calculating orbits very close to the pri- 
mary, and since these regions are in any case the least afl'ected 
by binarity, we consider a ring-like configuration for the initial 
disc: 0.6 < r < 1.1 (in units of rent)- 

For the binary, we consider a fixed mass ratio fi = 
0.25, corresponding to estimates of the mean mass ratio in 
binaries derived from extensive surveys (iKroupa et al.l Il990l: 
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Table 1. Set up for the disc-in-a-binary runs 



PARENT BODY RUN 




Number of test particles 


NpB = 2x10^ 


Initial radial extent^ 


0.6 < r< 1.1 


Initial eccentricity 


0.01 <eo< 0.05 


COLLISIONAL RUNS 




Average optical depth ^ 


(t> = 2 X 10-3 


Number of test particles 


Nnum = 2X10^ 


Size range ^ 


= 0.012 </3{s) </3(s,nin) = 2.5 


Size distribution at ^ = 


dN(s) oc s'^-^ds 



^ normalized to r^nY 
^ to derive focoii (Equ[T]) 
as parameterized by j3 oc \/s 



'Duquenno y and Mayoj 1 19911) . and explore 4 different orbital 
configurations: = 0, cb = 0.2, cb = 0.5 and cb = 0.75 (the 
values of gb are then automatically obtained from the require- 
ment that rcrit = 1). 

We consider an average vertical optical depth (r) = 2 x 10"^, 
typical of dense debris discs like /3 Pictoris. 

For the parent body runs, all particles start on circular or- 
bits (e = 0) in the binary's midplane (/ = 0). For the coUi- 
sional runs, particles are released from the parent bodies seeds 
following a size distribution in dN(s) oc s'^-^ds. Particle sizes 
are parameterized by the value of their JS parameter. Sizes are 
distributed continuously between jSmax = 0.5 and jSmin = 0.012. 
We take the minimum particle size to be fimax = 0.5, correspond- 
ing to the smallest grains on bound orbits. The contribution of 
smaller grains, below the blow-out size limit, is negligible for 
systems with optical depths in the (r) ~ 10"^ range (see Fig. 3 of 
iThebault et al.ll201Ql) . The number of sampled positions within 
a binary orbit is risav = 10. Convergence test runs have shown 
that this is the optimal value: higher values do not yield signif- 
icant changes in the final recombined density maps, while runs 
with lower risav lead to diverging results. Note that for the first 
binary orbit of the collisional run, the sample value is higher, 
^sav = 100, in order to accurately track the rapid initial blow-out 
of high-yS grains. All initial parameters are summarized in Tab[T] 

3.2. Results 

Depending on the considered configuration, the parent body sim- 
ulations are run for 200 (cb = case) to 500 (^g = 0.75) orbital 
periods of the binary in order to reach a steady state. This steady 
state is reached once all particles on unstable orbits have been 
removed and once all transient dynamical feat ures h ave disap- 
peared. As shown by Augereau & Papaloizoul (|2004|) . the most 
notable of these transient features are spiral arms, which de- 
velop because of the sudden introduction of a perturbing body, 
and disappear on the scale of a few hundred orbital periods. We 
then perform the collisional runs following the procedure pre- 
sented in Sec j2.2l The final disc profiles, after recombination of 
all collisional runs following the method described in Sec l2.3l 
are displayed in FigJS] 

A first important p oint is that these results confirm the main 
conclusion of Theba ult et al.l |2010) that circumprimary debris 
discs extend far outside the outer limit for orbital stability. For 
all four binary configurations the regions beyond rcrit = 1 are 
populated by small grains steadily produced by collisions in 
the parent body ring and placed on eccentric orbits by radia- 
tion pressure. The level of dustiness in these "forbidden" re- 
gions increases with cb because, for highly eccentric binaries. 



the perturbing companion is far from the main ring between two 
passages at periastron, so that small grains placed by radiation 
pressure in the r > rcrit region have a longer dynamical survival 
timescale before being removed. 

Regarding the spatial structure we find that, for low cb val- 
ues, precessing spiral features directed towards the companion's 
position develop from the main ring 0. These precessing spirals 
are the most pronounced for the = case, being clearly vis- 
ible for the whole binary orbit with no significant shape evo- 
lution, as is logically expected for this circular orbit case. As cb 
increases, these spiral structures become fainter. For the = 0.2 
case, they are only visible during half of the binary orbit (from 
periastron to apoastron), while they only briefly appear at perias- 
tron passages in the = 0.5 run. For the = 0.75 case no spi- 
ral is longer visible. This disappearance of the precessing spirals 
for high eB reflects the more general trend of the disc tending 
towards a time-invariant structure. In particular, no significant 
shape change is observed at diff'erent orbital phases of the bi- 
nary. The disc assumes a fixed and strongly asymmetric shape, 
extending much further out in the apoastron direction than in the 
periastron one (see radial profiles in Fig|4t and g). 

Another important result is the presence of strong inhomo- 
geneities in the main parent body ring itself. For low cb cases, 
they consist of a pronounced dip in the opposite direction of the 
companion's position (Fig|4]3), which precesses at the binary's 
angular velocity. In the moderately eccentric case (cb = 0.2), 
this dip is surrounded by two overdensitites, precessing at the 
same rate and symetrically positioned at ±70 - 90 with respect 
to the companion's antipolar position (Fig|4]l). The density con- 
trast around the density dip can reach up to ~ 60% (Fig Jit). For 
the cb = 0.2 case, in addition to the precessing features, a time- 
invariant structure appears, i.e., a dip in the binary orbit's peri- 
astron direction (Fig|4]l). As cb increases even more (cb > 0.5), 
this time-invariant asymmetry gets more pronounced, while the 
precessing structures progressively fade away (Fig|4f and h). 

3.3. Discussion 

The results of the previous section show that a companion star 
can significantly aff'ect the shape of a circumprimary disc, in- 
ducing pronounced radial and azimutal asymmetries. Basically, 
two diff'erent regimes can be identified depending on the binary's 
eccentricity. 

For low cb binaries, the disc's shape is time- varying and 
precessing with the binary's orbit. Strong asymmetries are ob- 
served, whose position is always related to that of the compan- 
ion on its orbit. Strong asymmetries are also observed for highly 
eccentric binaries, but they are in this case almost time-invariant 
and the disc has roughly the same shape regardless of the loca- 
tion of the companion on its orbit. In these high cb cases, the 
main asymmetries are oriented with respect to the binary orbit's 
geometry, not to the actual position of the companion on the or- 
bit. Between these 2 regimes, there is an intermediate state (see 
the = 0.5 runs) where the disc has a globally invariant asym- 
metric structure with transient features close to the companion's 
periastron passages. 

These fixed, precessing or transient features are caused by 
the combination of 3 distinct processes: 1) Collisional activity 
within the PB ring, steadily producing vast quantities of small 



3 These precessing spirals outside the parent body region should not 
be confused with the transient spirals that develop within the parent 
body ring and disappear once the initial conditions have been relaxed 
(see previous paragraph) 
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Fig. 3. Normalized vertical optical depth, at collisional and dynamical steady state, for 4 binary configurations; from top to bottom: 
cb = 0, cb = 0.2, cb = 0.5 and cb = 0.75. For each case, the disc profile is shown at periastron (left) and apoastron (right) passages 
of the companion star. The dotted circle is the theoretical outer radius for orbital stability rcru- The rectangle in the upper right corner 
shows the orbit of the binary as compared to the rcru radius (dotted circle) and the current position of the companion on its orbit. 
(An animated version of these results, displaying disc profiles for 10 diff'erent orbital positions of the companion, can be found at 
^http ://lesia.obspm.fr/perso/philippe-thebault/bindeb.html) . 
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Fig. 4. Steady-state surface density profiles for the same (from top to bottom) = 0, = 0.2, = 0.5 and et, = 0.75 cases as 
in FigJS] Left panels: radial profiles, along the binary semi-major axis direction, at periastron (full line) and apoastron (dotted line) 
passages of the companion star. Right panels: azimutal profiles, for a Ar = 0.1 wide annulus centered around r = 1 (in rcru units), at 
periastron (full line) and apoastron (dotted line) passages of the companion star. 
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grains that carry a large fraction of the geometrical cross sec- 
tion, 2) Radiation Pressure on the these small fragments, which 
places them on high-eccentricity orbits, 3) Gravitational pertur- 
bations by the companion star that will eventually, but not im- 
mediatly, remove all grains beyond rent- Equilibrium between 
these 3 effects results in a steady- state where a large fraction of 
the small, high-/? fragments produced in the PB ring remain in 
the dynamically forbidden regions long enough to cause the ap- 
pearance of pronounced structures, whose shape and evolution 
strongly depends on ^g. 

A careful examination of the variety of structures obtained 
for the different cases displayed in Fig|3] shows that the action 
of a companion star could be a potential explanation for some 
features that have been routinely observed in debris discs: spi- 
ral arms, azimutal inhomogeneties in ring-like structures, blobs, 
etc. It is worth stressing that all these structures are obtained at 
steady -state, and thus do not appear during a brief specific pe- 
riod in the system' s history, as would be the case for a fly-by 
with a passing star (lArdila et al.ll2QQ5l:lReche et al.ll2QQ9l) . 

Of course, this explanation should be taken with some cau- 
tion. A first reason is that other mechanisms have already been 
found to be able to produce asymmetries and inhomogeneities 
in discs, the m ost commonly invoked being the gravitational 
pull of a planet (lKrivovll2QlQl) . Another important point is that, 
contrary to the planet-perturbation scenario for which the pres- 
ence of an undetected planet can often be freely postulated 
given the observational constraints on the system, the presence 
or absence of a companion star is usually much better con- 
strained. Of all the 25 resolved debris discs (as of August 2011, 
see http://circumstellardisks.org/), only 3 are known to inhabit 
a confirmed binary system: HR4796A, HD181296 and f Ret 
(HD 14 1569 A is more a "transitional" disc than a bona fide de- 
bris di sc in the sense of the definition given by Lagrang e et al.l 
(120001) ). A detailed analysis of these systems exceeds the scope 
of the present numerical paper, as it would require exploring sev- 
eral additional free parameters, such as r or the mass ration jd. It 
will be undertaken in a forthcoming study. 

4. Limitations and Perspectives 

It should be pointed out that this code is only a first, albeit im- 
portant, step towards a fully integrated dynamics -h collisions 
model. Let us here list its main limitations and the possible im- 
provements that can be expected. 

- The most obvious limitation is that collisions are not mod- 
elled in a self-consistent way, but through an empirically de- 
fined collision destruction probability assigned to each par- 
ticle, and that no feedback of the collisions on the dynamics 
is included. These two problematic issues can probably not 
be handled by an N-body approach, for which the present 
code is maybe the upper limit of what can be achieved, at 
least for the case of high- velocity fragment-producing colli- 
sions. For this case of fragmenting impacts, a fully integrated 
and self-consistent dynamics-Fcollisions model is still out of 
reach today, and can probably only be achieved by the next 
generation of "hybrid" codes i n the spirit of the pioneering 
model explored by lGrigorieva et al.l (l2007h . 

- Our procedure is suited for systems where a level of symme- 
try has been reached, which is that the location of the parent 
bodies is symmetric with respect to the perturber's orbit, i.e., 
the global shape of the PB disc is identical between 2 pas- 
sages of the perturber at the same orbital position. Note that 
this symmetry is not required for each individual PB's orbit. 



which would be impossible because of the different period- 
icities between PBs and the perturber. What is required is a 
global symmetry for orbital streamlines, providing that each 
parent body's orbit is populated by a large number of bod- 
ies distributed uniformly in mean anomaly. This assumption 
of a phasing between global PB locations and pertuber orbit 
might not be valid for all perturbed disc configurations, in 
particular for systems strongly shaped by mean motion reso- 
nances. This case, usually occuring for discs interacting with 
embedded planets, will be explored in a forthcoming study. 
Note however that, for the present case of a circumprimary 
disc in a binary, this global symmetry is not assumed or pos- 
tulated but is a verified behaviour of the PB disc, which al- 
ways reaches this state after a transitory period of 100-300 
binary periods (see Sec l3.2b . 

- The equation governing collisional probability does not re- 
solve vertical structure in the disc because it is based on the 
vertically-integrated optical depth. This restrics the current 
version of the model to 2 dimensional problems, but those 
have been by far the most widely investigated by all past 
studies. A 3-D version of the collisional probability prescrip- 
tion is however fully implementable, at the cost of an in- 
crease in the size of the saved density maps at the end of the 
parent body runs. 

- Our model assumes that all collisions take place in the region 
of the parent body belt. This is not a major limitation for the 
present case of highly collisional massive debris discs, but 
the model will not perform well for fainter, drag-dominated 
discs where collisions can occur far interior to the belt. Our 
code can thus not include significant drag effects in its cur- 
rent version. 



5. Conclusions 

We have developed a new code to study the spatial structure of 
dynamically perturbed debris discs (be it by planets or compan- 
ion stars) taking into account the effect of collisions. It is espe- 
cially aimed at quantifying the competing effects between steady 
collisional production of small grains, placed by radiation pres- 
sure in dynamically unstable regions, and removal by gravita- 
tional perturbations. The principle of our numerical method is 
the following: we assume that the system has reached steady 
state and we perform a series of distinct N-body runs, each cor- 
responding to a different initial position of the perturber on its 
orbit, for which each particle has a collisional destruction prob- 
ability depending on its size and location. For each run, all par- 
ticle positions are recorded at regularly spaced time intervals. 
This database of "collisional" runs is then used to reconstruct 
the steady state profile of the disc at different orbital phases of 
the perturber. 

To illustrate the potential of this numerical procedure, we 
have applied it to the case of a circumprimary debris disc in 
a binary. We show that the complex interplay between colli- 
sional activity, radiation pressure and gravitational perturbations 
can create pronounced structures inside and outside the dynam- 
ical stability region. Depending on the binary's eccentricity, two 
regimes can be distinguished. For low ^g, the disc's structure is 
time varying, with spiral arms forming in the dynamically "for- 
bidden" region beyond rent and precessing at the binary's an- 
gular velocity. For high ^g, spiral structures disappear, the disc 
adopts a time invariant structure, extends far outside the stability 
region in the binary's apoastron direction and has a pronounced 
material depletion, inside rcrit, in the periastron direction. 
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Our model has the potential to be applied to many other con- 
figurations, in particular debris discs sculpted by planets, which 
will be the purpose of a future study. 
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